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Abstract-  The  presented  implementation  is  a  three-dimensional 
isotropic  monodomain  reaction-diffusion  model  with  a  realistic 
geometry  coupled  with  the  atrial  ion  model  proposed  by 
Nygren  et  al  The  partial  differential  equations  are  solved  by  a 
Garlekin  finite  element  method  in  space  and  a  forward  Euler 
approximation  in  time.  Simulations  yields  the  expected  results 
and  the  computational  performance  of  the  model  is  good 
considering  the  size  of  the  problem. 
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I.  INTRODUCTION 

A  fine  granularity  computer  model  would  be  useful  to 
further  elucidate  mechanisms  behind  atrial  arrhythmia.  This 
kind  of  information  is  hard,  if  not  impossible,  to  extract 
from  animal  models  and  in  vitro  experiments.  Knowing  the 
involved  ionic  currents,  new  pharmaceuticals  acting  on  the 
responsible  channels  might  be  developed.  The  possibility  of 
simulating  activation  patterns  may  be  useful  in  determining 
the  electrophysiological  effect  of  surgical  procedures. 
Maybe  pacing  strategies  can  be  developed  that  restore  the 
atrial  rhythm. 

Two  approaches  towards  simulation  of  the  cardiac 
electrical  activity  constitute  the  mainstream  of  proposed 
models;  cellular  automata  models  and  reaction-diffusion 
models  [1,2].  Reaction-diffusion  models  assume  that  the 
pseudo-syncytical  structure  of  the  cardiac  tissue  can  be 
regarded  as  a  single  volume  and  uses  a  diffusion  equation  to 
model  passive  spread  of  the  action  potential  and  an  ion 
current  model  represents  the  active  phase  (the  “reaction”). 

The  existing  ion  models  of  cardiac  myocytes  are  quite 
complex.  There  are  several  ionic  currents  and  cellular 
subspaces  involved  in  the  formation  of  action  potentials  and 
the  recovery  afterwards. 

Simulation  using  reaction-diffusion  models  require 
numerical  methods.  Traditionally  finite  difference  methods 
have  been  applied.  However,  Rogers  and  McCulloch  have 
presented  a  hybrid  finite  element  method  [3]. 

II.  Methodology 

The  implemented  model  is  an  isotropic  monodomain 
reaction-diffusion  model [2,4]  with  von  Neumann  boundary 
conditions  coupled  with  the  membrane  current  model  of 
Nygren  et  al  [ 5].  The  equation  to  be  solved  is  thus 


D  V  V  •  n  =  0  at  boundary 


where  V  is  the  transmembrane  voltage,  Iion  the  total 
sarcolemmal  ion  current  calculated  from  the  Nygren  model 
and  Ism  stimulation  pulses.  D  is  the  conductivity  tensor,  Cm 


the  transmembrane  capacitivity  and  /?  the  cell  surface  area  to 
volume  ratio,  n  is  a  vector  orthogonal  to  the  boundary. 

The  conductivity  tensor,  D,  is  in  the  anisotropic  case  a 
second  rank  tensor  but  reduces  to  an  invariant  in  an  isotropic 
model. 

A.  Solution  method 

To  generate  the  element  matrices  we  will  apply  the 
Garlekin  finite  element  method  to  the  reaction-diffusion 
equation. 

The  domain  (the  myocardium)  is  subdivided  into  a  set 
of  finite  cubic  elements  with  nodes  located  at  the  vertices.  In 
this  paper  we  use  a  Cartesian  global  coordinate  system  and 
element  coordinate  systems  aligned  with  the  global  system 

(Fig  1). 


Fig.  1 .  Definition  of  the  local  coordinate  system  for  element  k. 


The  value  of  V,  Iion  and  Ism  within  the  elements  is 
interpolated  from  the  value  of  the  corresponding  variable  at 
the  nodes.  A  set  of  linear  Lagrangian  interpolation  functions 
are  defined  in  the  element  coordinate  system: 

Vo  =<Po(^l)9ofe)<Pofe)  V4  =  <Podl)<Pofe>)<Pl(^3) 

Vi  -<Pi(^i)<Pofe-)<Pofe)  V5  =91(^1)90^2)91(^3) 

V2  =  <Po(^i)<Pi(^2)<Pofe)  Ve  =  <Po(^i)<Pi(^2)<Pifes) 
V3=<Pl(^l)<Pl(^2)<Po(^3)  V7=<Pl(^l)<Pl(^2)<Pl(^3) 
where 
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To  approximate  V  we  use  a  weighted  residual 
formulation  [6]  of  the  reaction  diffusion  equation,  here 
written  in  tensor  notation  (summation  convention  applies), 
with  the  interpolation  functions  as  weighting  functions.  We 
will  also  apply  the  product  rule  of  derivatives  followed  by 
the  divergence  theorem  to  the  left-hand  side,  the  surface 
integral  term  will  then  vanish  due  to  boundary  conditions. 
Finally,  we  rearrange  the  equations  and  move  nodal 
parameters  (since  they  are  independent  of  space  coordinates) 
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and  (3  and  Cm  (they  are  assumed  to  be  constant  over  the 
entire  domain)  out  of  the  integral: 

JJ|  r fk\aTfk\b  4i4243  = 

dt  k= 10  0  0 

M \\\DqPTt(k\a,r^bj(k\b,s^d^2d^- 
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where  denotes  a  mapping  from  index  a  in  element  k  to 
the  global  index  i.  Referring  to  Fig  1  we  note  that 

dxq  [o  p 

The  integrals  may  be  evaluated  symbolically  and 
assembled  into  the  global  matrices,  i.e.  map  the  local  node 
numbers  to  their  corresponding  global  node  numbers  and 
sum  the  contributions  from  each  element.  The  result  is  a 
system  of  ordinary  differential  equations  (ODEs): 
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one  term  of  the  sum  is  deduced  each  iteration.  This  is 
more  efficient  than  performing  the  transposition  and  the 
equivalence  is  easily  seen  in  the  formula  stated  above. 
Presentation 

The  visualizer  draws  a  three-dimensional  picture  using 
the  mesh.  The  nodal  potentials  are  then  color-mapped  to 
show  the  potentials  at  a  given  time.  By  selecting  a  point  a 
graph  of  the  potential  variation  over  time  is  drawn. 

C.  Simulation  parameters 

A  voxel  image  of  the  human  atria  was  provided  by  Rok 
Hren  at  the  University  of  Utah  (courtesy  of  Milan 
Horacek)[8].  The  final  mesh  consists  of  58235  nodes  and 
40681  elements. 

In  the  presented  simulations  the  ion  model  was  left 
unaltered,  leaving  only  the  capacitivity  Cm,  the  conductivity 
D  and  the  area  to  volume  ratio  /?  to  be  assigned.  In  this 
application,  due  to  the  fact  that  the  model  is  isotropic,  D  is  a 
scalar. 

Nygren  et  al.  define  Cm  for  the  whole  cell  to  be  0.050  nF 
and  Iion  is  defined  as  total  current  passing  into  the  cell;  no 
values  regarding  the  assumed  area  are  presented.  They  do, 
however,  assume  an  intracellular  volume  (Vo l,)  of  0.005884 
nL.  Using  these  values  we  may  easily  calculate  all  the 
necessary  parameters  except  D.  Since  the  presented  values 
for  the  intracellular  conductivity  varies  considerably  in  the 
literature  we  have  chosen  to  adjust  D  to  give  a  reasonable 
conduction  velocity;  to  simulate  normal  activation  D  was  set 
to  2.0  mS/mm.  This  is  within  the  range  of  presented  values. 

Stimulation  pulses  with  an  amplitude  of  1.5  nA  and  a 
duration  of  6  ms  was  added  to  the  ion  currents  of  the  nodes 
of  two  surfaces  (i.e.  to  eight  nodes)  close  to  each  other  to 
invoke  action  potentials. 

III.  Results 

A.  Computational  performance  and  memory  requirements 

The  routine  of  primary  interest,  the  one  that  calculates 
the  ion  currents  and  repeatedly  solves  the  equation  system, 
requires  about  12  hours  to  simulate  1  second  of  atrial 
activity. 

The  voxel  image  file  and  the  mesh  file  require  2  and  3 
MB,  respectively.  The  assembled  matrices  contains  1.5 
million  entries  each  and  the  factorized  8.5  million  entries. 
The  storage  required  by  these  matrices  adds  up  to  130  MB. 
The  state  file  contains  the  potential  and  the  28  ion  model 
state  variables  at  the  each  node;  this  requires  13  MB. 

The  size  of  the  simulation  file  (i.e.  the  node  potentials  at 
selected  times)  varies  due  to  the  run-length  encoding 
scheme;  as  a  rule  of  thumb,  200  MB  per  simulated  second. 

B.  Simulation  results 

Fig  3  shows  the  potential  distribution  from  a  “top  view” 
at  selected  times  after  a  stimulation  pulse  was  applied  in  the 
right  atrium  at  point  between  vena  cava  inferior  and 
superior.  Fig  3  also  shows  the  action  potential  at  a  selected 
point  in  the  model. 


Fig.  3.  Activation  pattern  following  stimulation  pulse.  This  figure  also 
shows  an  action  potential  at  the  marked  position. 


IV.  Discussion 

A.  The  model 

The  implemented  ion  model  is  the  one  proposed  by 
Nygren  et  al  This  is  the  most  recent  published  model  of  the 
atrial  myocyte  and  is  based  mainly  on  human  data.  In  the 
future  a  SA-model  will  be  used  for  the  SA-node.  This  model 
incorporates  the  ion  currents  and  concentrations  believed  to 
be  most  significant.  These  ion  currents  and  concentrations 
are  altered  in  ischemic  and  metabolic  conditions  and  by 
modifying  the  magnitude  and  timing,  of  these  currents  and 
concentrations  the  model  ought  to  be  able  to  reflect  many 
clinical  situations.  Hormones  and  neurotransmitters 
modifying  ion  channel  characteristica  may  also  be 
incorporated.  Structural  properties,  such  as  infarctions,  may 
be  incorporated  in  the  mesh  or  the  diffusion  tensor. 
Pharmaceutica  usually  either  acts  as  neurotransmitter 
agonists  or  antagonists  or  as  blockers  of  specific  ion 
channels.  Again,  by  modifying  ion  channel  properties,  the 
effect  of  these  pharmaceutica  may  be  simulated. 

Surgical  procedures,  such  as  maze  operations  and 
catether  ablations,  modifies  the  (electrical)  structure  of  the 
atria  by  introducing  scars  or  infarctions.  These  structural 
changes  may  be  reflected  as  changes  in  the  diffusion  tensor. 
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Electrical  therapies  may  be  modeled  by  introducing 
stimulation  pulses  at  selected  points  in  the  model.  Input 
from  sensing  electrodes  is  available  since  the  potential 
distribution  is  known.  Pacing  algorithms  are  easily 
interfaced  to  existing  code. 

Due  to  lack  of  data,  the  implemented  model  is  isotropic. 
Though,  extending  it  to  an  ansisotropic  model  would  only 
require  an  estimate  of  the  diffusion  tensor  and  the  addition 
of  a  quadrature  algorithm 

B.  Solution  method  and  implementation 

Finite  element  methods  may  easily  handle  von 
Neumann  conditions  when  Garlekin  weighting  functions  are 
used;  the  complexity  of  the  problem  is  actually  reduced.  The 
explicit  method  used  to  approximate  the  time  derivatives  of 
potentials  as  well  as  in  the  differential  equations  of  the  ion 
model  poses  a  risk  of  instability.  The  required  time  step 
might  be  smaller  than  what  would  be  required  by  an  implicit 
method,  but  we  believe  that  the  speed  of  the  explicit  method 
most  likely  outweigh  the  reduction  of  time  step  possible 
with  an  implicit  method. 

The  design  criteria  may  summarized  as:  accuracy , 
speed ,  extensibility  and  portability.  These  apply  in  variable 
degree  to  different  steps  in  the  solution  process.  Accuracy  is 
mostly  governed  by  the  choice  of  solution  method  and  ion 
model;  at  the  implementation  level  there  are  three  factors 
that  needs  to  be  considered:  floating  point  precision,  order 
of  expression  evaluation  and  the  magnitude  of  the  time-step 
taken  during  the  solution  of  the  ODEs.  Speed  is  a  main 
concern  when  the  equation  system  is  to  be  solved  due  to 
large  number  of  iterations,  while  longer  run-times  are 
acceptable  in  the  preparation  of  the  matrices.  Also,  if  the 
run-time  is  rendered  unacceptable,  emigration  to  a  more 
powerful  platform  would  be  of  interest.  Thus,  portability  is 
of  importance. 

C.  Computational  performance  and  memory  requirements 

The  present  computer  (Pentium  III  800MHz)  requires 
about  2  seconds  for  each  timestep  due  to  the  large  number  of 
equations,  and  the  expensive  exponential  operations 
required  by  the  nonlinear  right-hand  sides  of  the  ion  model. 
This  accumulates  to  a  barely  acceptable  simulation  time  of 
roughly  12  hours  for  each  simulated  second. 

There  are  today  two  steps  in  the  solution  process  that 
are  outweighing  the  others;  the  two  triangular  solve  steps 
(forward  elimination  and  backward  substitution)  and  the  ion 
model  operations.  The  factorization  step  generates  a 
triangular  matrix  containing  about  8.5  million  entries  with 
the  permutation  algorithm  used  today.  The  current 
permutation  algorithm  is  balanced  between  the  number  of 
fill-ins  generated  and  the  time  used  to  find  the  permutation 
considering  a  system  that  is  supposed  to  be  solved  once. 
This  is  far  from  optimal  in  this  application  since  the 
triangular  solve  is  iterated  fifty  thousand  times  for  each 
simulated  second  while  the  permutation  is  only  sought  once. 

We  also  plan  to  develop  a  parallelized  version  of  the 
algorithm. 


D.  Simulation  results 

Simulations  of  a  single  wavefront  performed  with  this 
model  show  a  reasonable  agreement  with  available  mapping 
studies [9, 10]  regarding  average  conduction  velocity.  The 
refractory  period  is  also  in  line  with  published  data.  The 
model  behaves  as  expected  when  two  activation  fronts 
collide. 

V.  Conclusions 

The  presented  implementation  is  an  three-dimensional 
isotropic  monodomain  reaction-diffusion  model  with  a 
realistic  geometry  coupled  with  the  atrial  ion  model 
proposed  by  Nygren  et  al.  Thus,  it  allows  conduction 
velocity  and  refractory  periods  to  adjusted.  Normal 
anisotropy  is  not  included  in  the  current  model  due  to  the 
absence  of  data  but  the  modification  is  rather  easy. 

The  computational  performance  of  the  model  is  good 
considering  the  size  of  the  problem.  However,  run-times 
needs  to  further  reduced  for  the  model  to  reach  its  full 
potential. 
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